Cone-beam computerized tomography (CT) reconstructs the density function of a three-dimensional object from a set of cone-beam projections. Such a system uses an area detector to receive rays emitted from an X-ray point source and attenuated by partial absorption in the object that they pass through. As in traditional (i.e., planar) CT, the source 10 and the detector 14 are placed on opposite sides of the object 12 being scanned (see FIG. 1). Rays contributing to an image on the detector surface form a cone with the X-ray source 10 at the apex. From the X-ray radiance value recorded at a point on the area detector 14, one can compute the integral of attenuation along the ray from the X-ray source 10 to the given point on the detector 14.
As the source-detector 10/14 pair undergoes a simultaneous rotation and translation around the object 12, a plurality of two-dimensional cone-beam images projected from various source positions can be acquired and used to reconstruct the distribution of absorption inside the three-dimensional object 12.
Compared to the traditional slice-at-a-time tomographic machine, the cone-beam CT offers faster scans, higher patient throughput, significant reduction in X-ray dosage, and isotropic resolution. It has a great potential to be applied to a wide range of medical and industrial applications.
Radon's 1917 inversion formula (Johann Radon “Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten,” Ber. Verh. Sächs. Akad. Wiss. Leipzig. Math. Nat. Kl., Vol. 69, pp. 262-277, 1917) plays an important role in understanding the cone-beam reconstruction problem. The building blocks of the three-dimensional Radon inversion formula are planar integrals. We can write a plane in R3 asLl,β:={x ∈Rn|x ·β=l, l≧0,βεSn−1},  (1)where β is the unit normal of the plane and 1 is the perpendicular distance of the plane from the origin. The Radon transform of a function f on R3 is defined as the set of integrals of ƒ over all the planes in R3 which can be expressed as a function of two parameters (l and β):                                           f            ⁡                          (              x              )                                =                                                    -                                  1                                      8                    ⁢                                          π                      2                                                                                  ⁢                                                ∫                                      S                    2                                                  ⁢                                                                                                    ∂                        2                                            ⁢                      R                                        ⁢                                                                                   ⁢                                          f                      ⁡                                              (                                                  l                          ,                          β                                                )                                                                                                  ∂                                          l                      2                                                                                            ⁢                          |                              l                -                                  x                  ·                  β                                                      ⁢                          ⅆ              β                                      ,                  x          ∈          Ω                ,                            (        3        )            
The Radon formula is given by:                               Rf          ⁡                      (                          l              ,              β                        )                          :=                              ∫                          x              ∈                              {                                                      x                    |                                          x                      ·                      β                                                        =                  l                                }                                              ⁢                                    f              ⁡                              (                x                )                                      ⁢                                          ⅆ                x                            .                                                          (        2        )            in which, S2 denotes the two-dimensional unit sphere in R3 and Ω denotes the support of ƒ. The integral in Eqn. (3) over S2 is the backprojection operator; it integrates over all the planes passing through x. The integration sphere is therefore called the backprojection sphere with its center at x, denoted by Sx2 (x is considered as an index). It is clear that the points on Sx2 represent the unit normals of all the planes through x.
To recover the function value at point x, R″ƒ(l,β) is obtained on all or almost all planes passing through x. In cone-beam reconstruction, however, planar integrals are not available from the cone-beam data because rays diverge from the point source inside each projection. Hence, the Radon formula (Eqn. (3)) is not immediately employed.
The first cone-beam inversion formula for real-valued functions is given by Tuy in 1983; this formula was a Fourier-based method (Heang K. Tuy “An Inversion Formula for Cone-Beam Reconstruction,” SIAM J. Appl. Math, Vol. 43, 1983, pp. 546-552). Smith's paper in 1985 established connections between the cone-beam data and the second-order radial derivative of the Radon transform, R″ ƒ(Bruce D. Smith “Image Reconstruction from Cone-Beam Projections: Necessary and Sufficient Conditions and Reconstruction Methods,” IEEE Trans. Med. Imag., Vol. 4, 1985, pp. 14-25). The most important contribution in these early derivations is a clear understanding of the data sufficiency condition for an exact reconstruction, that is, almost all planes passing by the support of the object shall intersect with the source orbit.
The next significant breakthrough came with the discovery of the Fundamental Relation by Grangeat (Pierre Grangeat “Mathematical Framework of Cone-Beam 3D Reconstruction via the First Derivative of the Radon Transform,” Mathematical methods in tomography, Lecture notes in mathematics 1497, 1991, pp.66-97. The Fundamental Relation relates the cone-beam data on a slice of fan-beams inside each cone-beam projection to the first-order radial derivative of the Radon transform, R′ƒ. R′ƒ serves as an implicit link between the cone-beam data and R″ ƒ. The second-order radial derivative of the Radon transform is needed in order to use Eqn. (3); R″ ƒis then backprojected to recover ƒ.
Though substantial progress has been made during the last two decades, the solutions for exact cone-beam reconstruction are still not filly satisfying. In many of the reconstruction methods that have been developed, the backprojection-differentiation operation inherited from the Radon formula appears ad hoc and is the most time-consuming step in the reconstruction.
The well-known filtered backprojection (FBP) cone-beam reconstruction technique, which is widely used in industry, is given by Feldkamp (FDK) et al. for circular source orbits (L. A. Feldkamp, L. C. David and J. W. Kress “Practical Cone-Beam Algorithm,” J. Opt. Soc. Am. A., Vol. 1, No. 6, 1984, pp. 612-619). In such a case, data from cone-beams with narrow angles is treated in an approximate way using extensions of two-dimensional fan-beam methods. The FDK algorithm is easy to implement; however, it only provides reasonably good reconstruction near the mid-plane and cannot be used for wide cone angles. Hence, alternative reconstruction methods and the embodying imaging apparatus are still being sought, particularly for the large-detector cone-beam system since it has become a reality.
In designing a dedicated cone-beam imaging system, finding a proper source orbit is a challenge. The selection of a good source orbit not only depends on the dimension of the object under investigation, but also depends on the geometric measurements such as the allowed source-to-object and detector-to-object separation. An important condition for accurate reconstruction is the data sufficiency condition. Another desirable feature of the selected source orbit is symmetry.
Among various source orbits that have been proposed, sinusoidal trajectory and helical trajectory meet both conditions. Though advantageous in their sampling performance, reconstruction procedure using these two scan paths have yet to achieve the desired efficiency. The principal difficulty encountered in the reconstruction is caused by the sophisticated mapping from the local projection geometry to the Radon space geometry characteristic to many non-planar source orbits.
Other approaches use two orthogonal planar trajectories such as circle-plus-circle, circle-plus-line and circle-plus-arc to fulfill the data sufficiency condition. Although the hybrid methods combine cone-beam data from two simpler scanning processes, they have two major disadvantages. First, the discontinuity in the mechanical movement makes them less attractive in practice. Second, sampling in the Radon space where the backprojection takes place is not balanced under these hybrid scanning geometries; this limits the reconstruction accuracy.